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The key feature for the successful implementation of the surrogate data test for nonlinearity on 
a scalar time series is the generation of surrogate data that represent exactly the null hypothesis 
(statically transformed normal stochastic process), i.e. they possess the sample autocorrelation and 
amplitude distribution of the given data. A new conceptual approach and algorithm for the gen- 
eration of surrogate data is proposed, called statically transformed autoregressive process (STAP). 
It identifies a normal autoregressive process and a monotonic static transform, so that the trans- 
formed realisations of this process fulfill exactly both conditions and do not suffer from bias in 
autocorrelation as the surrogate data generated by other algorithms. The appropriateness of STAP 
is demonstrated with simulated and real world data. 
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The surrogate data test for nonlinearity has been 
widely used in real applications in order to establish sta- 
tistically the existence of nonlinear dynamics and justify 
the use of nonlinear tools in the analysis of time series 
10, ||] . The most general null hypothesis Hq of this test so 
far is that the examined time series x = [xi, . . . , Xn]' is 
a realisation of a normal (linear) process s= [si, . . . , s„]' 
undergoing a possibly nonlinear static transform /i, i.e. 
Xi = h[si), i = 1, . . . ,n. To test the null hypothesis, an 
estimate q from a nonlinear method applied to the orig- 
inal data, goj is compared to the estimates, 91, . . . ,9a/ 
on M surrogate time series representing the null hypoth- 
esis A properly designed surrogate time series z 
should possess the same autocorrelation as the original 
data, rz{T) = r^ir) for a range of lags r, and the same 
amplitude distribution, Fz{zi) = Fx{xi) {Fx{xi) is the 
cumulative density function (cdf) of x^), and be other- 
wise random. 

It has been reported that erroneous results are likely to 
occur mainly due to the insufficiency of the algorithms to 
generate surrogate data that preserve the original linear 
correlations The prominent algorithm of amplitude 
adjusted Fourier transform (AAFT) Q, used in most 
applications so far, is built based on the assumption of 
monotonicity of h. When h is not monotonic, the AAFT 
algorithm is found to favour the rejection of Hq due to 
the mismatch of the original linear correlations 0]. The 
iterated AAFT (lAAFT) algorithm improves the match 
of the autocorrelation of AAFT ||^ , but with about the 
same accuracy for all the surrogates, so that the small 
variance, combined with the small bias, may be another 
cause for false rejections in some cases 0. Another algo- 
rithm making use of simulation annealing seems to per- 
form similarly to lAAFT Recently, a correction of 
AAFT (CAAFT) that results in unbiased match of the 
linear correlations was proposed in |l0| . 

In this paper, we develop further the somehow pro- 



found rationale of the CAAFT algorithm and formulate 
a new conceptual approach for the generation of surro- 
gate data consistent with Ho, which we then solve an- 
alytically. The new algorithm, called STAP, generates 
the surrogate data as realisations of a suitable statically 
transformed autoregressive process (STAP), i.e. the pro- 
cess under Hq is designed as a static transform of a suit- 
able normal process. 

The main idea behind the STAP algorithm is that for 
any stationary process X = [Xi, X2, ...]', for which we 
have finite scalar measurements x, there is a scalar lin- 
ear stochastic process Z with the same autocorrelation 
px and marginal cdf $x as for the observed process 
{Pzir) = Px{t) and ^z{zi) = ^x{xi)), i.e. Z is a scalar 
"linear copy" of the observed X. The objective is to de- 
rive Z through a static monotonic transform g on a scalar 
normal process U with a proper autocorrelation i.e. 
Zi ~ g{Ui). Thus g and pu have to be properly selected, 
so that Z has the desired properties. In practice, the sur- 
rogate data set is cL finite realiscition z — [zi, ■ . . , Zn Y of 
the process Z and g and pu are estimated based solely 
on X. Note that g and U are in general different from 
h and S, respectively, of Hq (they are the same if h is 
monotonic Thus with this approach Hq can be for- 

mulated more generally that the time series is generated 
by a linear stochastic process. 

Let $0 be the marginal cdf of a standard normal pro- 
cess U. A suitable choice for g, so that ^z{zi) — ^x{xi), 
is defined as lOl 



Z,=gm = ^x\MU^)), 



(1) 



* Electronic add ress : 
auth . gr/ dkugiu] 



lkugiu@gen.auth.gr; URL: tittp: //users 



where g is monotonic by construction. Assuming that 
^xixi) is continuous and strictly increasing and that 
~1 < Pu < 1, which are both true for all practical pur- 
poses, there is a function (j) dep ending on g, such that 
Pz = 4>{pu) for any lag r ||l2[ |l3|. If g has an ana- 
lytic form, then it may be possible to find an analytic 
expression for cj) as well. In that case, given that px is 
known and by setting pz = px, one can invert 4> to find 
Pu ~ 4>~^ipx), if 4>~^ exists. 

In general, the function g, as defined in eq.(p]), does not 
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have an analytic form because is not known analyti- 
cally, but it can be approximated by an analytic function, 
e.g. a polynomial Pm of degree m. 



g{Ui) ~ PmiUi) = Qq + 



(2) 



Low degree polynomials have been used to approximate 
such transforms |l^ . Then using the definition for the 
autocorrelation, the approximate expression for (j) reads 



Em v^m / \ 1 



(3) 



where an arbitrary lag r is implied as argument for the 
autocorrelations, Hs is the sth central moment of Ui being 
P-2k+i = 0, fj,2k = 1 ■ 3 ■ • • (2fc - 1), fc > 0, and ^s,t is the 
sth-tth central joint moment of the bivariate standard 
normal distribution of {Ui, Ui-r), defined as follows |^ 



(2p, 



l2j 



(2A:)!(20!y- 



{2k + l)\{2l + l)\^ 
/i2fc+i,2;+i - ^^^TiVi 

3=0 

fiij — if k + I = odd. 
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2j + l 



(fc-j)!G-j)!(2j + l) 



where v = min(fc, Z). By substituting the expression 
for the moments in eq.(|^), the expression for cj) can be 
brought to a polynomial form of the same order m 



Px = (t>{pu) 



(4) 



where the vector of coefficients c = [ci, . . . ,Cm]' is ex- 
pressed only in terms of a = [ai, . . . (the expres- 
sions are rather involved and therefore not presented 
here). Simpler expressions can be derived using the 
Tchebycheff-Hermite polynomials [ ^ . Thus an analytic 
expression for pu is possible if eq.(^can be solved with 
respect to pu- From our simulations, we conjecture that 
if g is monotonic then (j) is also monotonic in [—1,1]. Then 
f/i"^ exists and a unique solution for pu can be found from 
eq.(Q). The proper standard normal process U is com- 
pletely defined by pu and applying the transform g of 
eq.(|^) to the components of U, the "linear copy" Z of 
the given process X is constructed. 

Note that the solution for pu is given analytically from 
the polynomial approximation of g and it requires only 
the knowledge of the coefhcients a of the polynomial and 
the autocorrelation px- 

In practice, we operate with a single time series x 
rather than a process X and with the sample estimates 
Fx and for <i>x and px , respectively. The steps of the 
algorithm are as follows: 

1. Estimate the vector of coefficients a = [ai, . . . , am]' 
of the polynomial Pm from the graph of Xi — 
F~^{Fo{wi)), i.e. the graph of x vs w after their 
ranks are matched, where w = [wi,... is 
standard white normal noise. 



2. Compute c 
eqs.(|^). 



[ci, . . . , Cm]' for the given a from 



3. Find from eq.(|j) for the given c and Tx using 
the sample estimates r instead of p. The common 
practice is that the solution exists and it is unique. 
If this is not the case, repeat the steps 1-3 for a 
new w until a unique solution is obtained. 

4. Generate a realisation u of a standard normal pro- 
cess with autocorrelation r„. We choose to do 
this simply by means of an autoregressive model 
of some order p, AR(p). The parameters b = 
[bo,bi,... jbpY of AR(p) are found from r„ us- 
ing the normal equations solved effectively by the 
Levinson algorithm |jl^. The AR{p) model is run 
to generate u 



Ui+i = ^0 + ^ bj 



-(i-i) 



N(0,1). 



5. Transform u to z by reordering x to match the rank 
order of u, i.e. Zi — F^^(Fo(ui)). 

Note that u possesses the sample normal marginal cdf 
-fo and the proper r„, so that z possesses F^ — Fx, 
Tz — Tx, and is otherwise random, as desired. In prac- 
tice however, the equality = Tx is not exact and 
may vary substantially around Tx- Two possible reasons 
for this are the insufficient approximation of g in step 1 
and the inevitable variation of the sample autocorrela- 
tion of the generated u in step 4, which decreases with 
the increase of data size. The former is due to the limited 
power of polynomials in approximating monotonic func- 
tions and this shortcoming causes also occasional repeti- 
tions of the first steps of the algorithm as stated in step 
3 . The latter constitutes an inherent property of the 
so-called "typical realisation" approach (i.e. a model is 
used to generate the surrogate data) and cannot be con- 
trolled. However, less variation in the autocorrelation 
is achieved when the AR(p) model is optimised making 
the following steps, in the same way as for the CAAFT 
algorithm ]lO| |: 

1 . Apply the algorithm presented above K times and 
get 2},. . . ,2,^ surrogate time series. 



2. Compute 
to r„ M. 



, r.K and find the one, r.i closest 



3. Use the parameters b of the ^-repetition to generate 
the M surrogate data (steps 4-5 of the algorithm 
above) . 

The K repetitions above as well as the occasional rep- 
etitions of steps 1-3 of the first part of the algorithm 
may slow down the algorithm if the time series is long, 
but they have no impact on the principal function of 
the algorithm. Simply, some realisations of white noise 
w are discarded in the search of the parameters b of 
the most suitable AR model that generates the surrogate 
data (through the g transform). 

The free parameters of the STAP algorithm are the 
degree m of the polynomial approximation of g, the or- 
der p of the AR model, the number K of repetitions for 
the optimisation of AR(p) and the maximum lag Tmax, 
used to compare r^i, . . . ,rzK to Tx- Usually, a small m 
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(m < 10) is sufficient. For p, there is no optimal range of 
values but it may vary with the shape of r^, e.g. a slowly 
decaying may be better modelled by a larger p. In all 
our simulations, we set K — M — 40 and Tmax = P- 

The proper performance and superiority of STAP over 
AAFT and lAAFT was confirmed from simulations on 
different toy models. CAAFT was found to perform very 
similarly to STAP. We show in Fig. |^ comparative re- 
sults for AAFT, lAAFT and STAP for three represen- 
tative synthetic systems: the cube of an AR(1) process 
(sj+i — 0.3 -I- 0.8si -I- Ci.ei ~ N(0, 1), Xi = s?), the square 
of the same AR(1) [xi = s^), both being consistent with 
Ho, and the x variable of the Rossler system |19 , not 
consistent with Hq. For each system, we generate 100 
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FIG. 1: Percentages of rejections of Ho using as discriminating 
statistics the fit of Volterra polynomials from 100 realisations 
for each of the three cases in the three panels, as indicated. 
Three types of surrogates are used in each test as shown in 
the legend (for STAP m = 5, p = 5). The vertical gray line 
distinguishes the linear from the nonlinear statistics. 

time series of 2048 samples each and for each realisation 
we generate M = 40 surrogate data of each type. 

As discriminating statistics g* we choose the correla- 
tion coefficient (CC) of the fit with the series of Volterra 
polynomials of degree 2 and order u = 5. The polyno- 
mials for the first i terms, where z = 1, . . . , 5, are linear 
and for terms i = 6, . . . ,20 are nonlinear (see also [p^). 
To quantify the discrimination we use the significance 



S 
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for each polynomial of i terms, where is 



the statistic on the original data, (q*) and are the 
mean and standard deviation of the statistic on the 
M surrogate data. The null hypothesis Hq is formally 
rejected at the 0.05 significance level when S > 1.96, un- 
der the assumption of normality for the statistic q, which 
turns out to hold in general. The percentages of rejec- 
tions for each of the three systems are shown in Fig. |^. 
Very similar results were found when deciding the rejec- 
tion from the rank ordering of 9ii ■ • ■ j Qm- 

For the linear statistics, STAP gives consistently and 
correctly no rejections, i.e. unbiased match of the lin- 
ear correlations, whereas AAFT and lAAFT give very 
large percentages of rejections for all but the first case 
where AAFT gives about 5% rejections, as expected. For 
lAAFT, the rejections occur because ct^ is very small (10 
to 20 times smaller than for AAFT and STAP), though 
the bias Qq — (g*) is smaller than for STAP. 
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FIG. 2: (a) The correlation coefficient of the fit with Volterra 
polynomials on the original normal EEG data (black line) 
and 40 AAFT, lAAFT and STAP surrogates (gray lines) in 
the three panels, as indicated (for STAP m = 10, p = 10). 
(b) The same as in (a) but for the epileptic EEG. (c) The 
significance for the fits in (a) (upper panel) and in (b) (lower 
panel) . The vertical lines in the plots distinguish the linear 
from the nonlinear polynomials. 



For the nonlinear statistics, the feature of the linear 
statistics persists for the two first systems (consistent 
with Ho) and for all three algorithms, i.e. correctly no re- 
jections with STAP and erronous rejections with AAFT 
and lAAFT. We cannot explain why the polynomial fit 
for the lAAFT surrogates improves with the addition of 
nonlinear terms (see Fig. |l|). For the nonlinear system, 
STAP properly converges with the addition of few non- 
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linear terms to the correct 100% rejection level, which 
AAFT and lAAFT possessed already with linear statis- 
tics. 

Next, we verify the three algorithms on two human 
EEG data sets, one recorded many hours before an 
epileptic seizure accounting for normal brain activity and 
another recorded during an epileptic seizure. The epilep- 
tic EEG seems to exhibit a pattern of oscillations indicat- 
ing a more deterministic character than the non-epileptic 
EEG, and this constitutes a well established result in 
physiology. This is demonstrated also with the fit of 
Volterra polynomials in Fig. ^. The fit improves with 
the inclusion of the first nonlinear terms for the epileptic 
EEG but not for the non-epileptic EEG. These findings 
are confirmed statistically by the test with STAP surro- 
gate data while the results of the test with AAFT and 
lAAFT are more or less confusing. 

In particular, the Hq on the normal EEG is erroneously 
rejected with AAFT because the difference in the fit be- 
tween original and surrogate data is about the same for 
the linear and nonlinear polynomials (see Fig. ^ and c). 
The same test result is obtained with lAAFT for large 
nonlinear polynomials, whereas again there is significant 
difference in the linear fits between original and lAAFT 
surrogates (not easily discernible as both bias and vari- 
ance are very small). Using STAP surrogates, the Ho is 
not rejected for both linear and nonlinear fits. 

For the epileptic EEG, there is again a clear differ- 
ence in the linear fit between original data and AAFT 
surrogates and a smaller but equally significant differ- 
ence between original data and lAAFT surrogates (see 
Fig. ^ and c). The significance S for both AAFT and 
lAAFT increases with the addition of the first couple of 
nonlinear terms, much more for lAAFT due to the small 



variance of CC. However, the deviation in the linear fit 
does not support reliable rejection of Hq. On the other 
hand, using STAP surrogates, the Hq is properly rejected 
only for the nonlinear statistics and with high confidence 
(5 < 2 for the linear fit and 5 ~ 5 for the nonlinear fit). 

In general, the test with STAP surrogates tends to be 
more conservative, i.e. small discriminations are found 
less significant, as the data size decreases. For example, 
the test on 296 sunspot samples (for which a small leap 
of the polynomial fit with the addition of nonlinear terms 
was observed) gave rejection of Hq for AAFT and lAAFT 
but not for STAP (not shown here, see also jl^). How- 
ever, this should not be considered as a drawback of the 
STAP algorithm, as one expects that the power of the 
test reduces with the decrease of data size. 

A new algorithm that generates surrogate data for the 
test for nonlinearity has been presented, called statically 
transformed autoregressive process (STAP). The key fea- 
ture of STAP is that it finds analytically the autocorre- 
lation of an appropriate underlying normal process for 
the test. This is the main difference of the STAP al- 
gorithm from the corrected AAFT (CAAFT) algorithm, 
where the autocorrelation is estimated numerically. Both 
CAAFT and STAP algorithms do not suffer from the se- 
vere drawback of the AAFT algorithm, i.e. bias in the 
match of the original autocorrelation. The AAFT algo- 
rithm is essentially impractical for real applications be- 
cause it favours the rejection of Hq as a result of the 
bias in the autocorrelation. From the numerical simula- 
tions, it turns out that the lAAFT algortihm may also 
give small bias in the linear correlations, favouring also 
the rejection of Hq. On the other hand, the STAP algo- 
rithm performs properly and gives reliable rejections of 
Hq, only whenever this appears to be the case. 
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